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Abstract 

A generalization of the Vandermonde matrices which arise when the power ba- 
sis is replaced by the Said-Ball basis is considered. When the nodes are inside the 
interval (0, 1), then those matrices are strictly totally positive. An algorithm for 
computing the bidiagonal decomposition of those Said-Ball- Vandermonde matrices 
is presented, which allows to use known algorithms for totally positive matrices rep- 
resented by their bidiagonal decomposition. The algorithm is shown to be fast and 
to guarantee high relative accuracy. Some numerical experiments which illustrate 
the good behaviour of the algorithm are included. 
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1 Introduction 



Numerical computing with structured totally nonnegative matrices is a classi- 
cal subject in the field of numerical linear algebra which has recently received 
a renewed attention, as can be seen in the recent survey paper [8] , where sev- 
eral different classes of structured matrices are considered, among them totally 
positive matrices. 

Classically, a matrix is said to be totally positive if all its minors are nonnega- 
tive [11]. Consequently, the matrices with that property are also called totally 
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nonnegative matrices [TT], and this term is becoming more used in recent lit- 
erature. 

The fact that a nonsingular totally nonnegative (TN) matrix can be decom- 
posed as a product of nonnegative bidiagonal factors was used by Koev [2T][22] 
to develop several accurate algorithms for the general class of TN matrices. A 
detailed survey of several results related to TN matrices, including the bidi- 
agonal factorization, has been presented in 

Nevertheless, it must be stressed that the algorithms of Koev pT][22] start from 
the bidiagonal decomposition of a TN matrix A, which is stored in a matrix 
which is denoted there as BV{A), and that such a decomposition needs to be 
computed for each particular class of TN matrices being considered. Using the 
words of the section devoted to conclusions and open problems in |22] : 

The caveat in our algorithms is that every TN matrix must be represented by 
its bidiagonal decomposition. While every TN matrix intrinsically possesses 
such a decomposition, and for many classes of structured matrices this decom- 
position is very easy to obtain accurately , there are important TN matrices for 
which we know of no accurate and efficient way to compute their bidiagonal 
decompositions. 

Examples of totally nonnegative matrices for which there are accurate and ef- 
ficient algorithms for computing BT>{A) are Vandermonde [Hll6fl7fl8j . Cauchy 
[S] , Cauchy- Vandermonde , generalized Vandermonde [TU] and Bernstein- 

Vandermonde matrices j25] . 

On the other hand, it is not always recognized that while Neville elimination 
[TT|12fl3P^ is a key theoretical tool for the analysis of that bidiagonal de- 
composition, it generally fails to provide an accurate algorithm for computing 
BV{A). This fact is explicitly noted in [20], where the author indicates that 
the function TNBD is the only function in the package TNTool that does not 
guarantee high relative accuracy. 

Consequently, the accurate (and, if possible, fast) computation of BV{A) is a 
previous task to be performed before applying Koev's algorithms to a given 
class of TN matrices. The importance of those algorithms was very acknowl- 
edged in [27], while relevant previous results were presented in [9]. 

In this work we are extending to the class of Said-Ball- Vandermonde matrices 
the work we have recently carried out for the class of Bernstein- Vandermonde 
matrices [23|. A crucial fact for obtaining high relative accuracy in our algo- 
rithm is that it satisfies what is called in [S] the NIC (no inaccurate cancella- 
tion) condition: 

NIC: The algorithm only multiplies, divides, adds (resp., substracts) real 
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numbers with like (resp., differing) signs, and otherwise only adds or substracts 
input data. 

The Said-Ball basis is a generalization of the Ball basis pl2]l3] . a well-known 
basis for cubic polynomials on a finite interval which is useful in the field of 
Computer-Aided Design. The Said-Ball basis was introduced for odd degree 
polynomials by Said in [26], and then its definition for polynomials of even 
degree was suggested in [TH] . Its properties in connection with total positivity 
and shape preservation were studied by Goodman and Said for odd degree 
polynomials [15], and recently by Delgado an Pena in [7], where it was estab- 
lished that the Said-Ball basis is a normalized totally positive (NTP) basis for 
every value of the polynomial degree. 

The rest of the paper is organized as follows. Some basic results on Neville 
elimination and total positivity are recalled in Section 2. In Section 3 the bidi- 
agonal decomposition of a Said-Ball- Vandermonde matrix and of its inverse 
are presented. The algorithm for computing these bidiagonal factorizations is 
introduced in Section 4. In Section 5 the problems of linear system solving and 
eigenvalue computation for a Said-Ball- Vandermonde matrix are considered. 
Finally, Section 6 is devoted to illustrate the accuracy of our algorithms by 
means of some numerical experiments. 



2 Basic results on Neville elimination and total positivity 



In this section we will briefly recall some basic results on Neville elimination 
and total positivity which we will apply in Section 3. Our notation follows the 
notation used in [T2p3] . Given k,nEJSi{l<k<n), Qk,n will denote the set 
of all increasing sequences of k positive integers less than or equal to n. 

Let y4 be a real square matrix of order n. For k < n, m < n, and for any 
« £ Qk,n and P G Qm,n, we will denote by ^[q;|/3] the submatrix k x m of A 
containing the rows numbered by a and the columns numbered by (3. 

The fundamental tool for obtaining the theoretical results applied in this paper 
is the Neville elimination (see [T^fT^ ). a procedure that makes zeros in a 
matrix adding to a given row an appropriate multiple of the previous one. For 
a nonsingular matrix A = {cii,j)i<i,j<n, it consists on n — 1 steps resulting in 
a sequence of matrices A := A\ A2 — >■ ... An, where At = i<ij<n 
has zeros below its main diagonal in the t — 1 first columns. The matrix A^^i 
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is obtained from At {t = 1, . . . ,n) hj using the following formula: 



a, 



''id 



if i<t 
, if i > t + 1 and j > t + 1 



(2.1) 



otherwise. 



In this process the element 



'id 



is called pivot of the Neville elimination of A. The process would break 
down if any of the pivots pij {j < i < n) is zero. In that case we can move 
the corresponding rows to the bottom and proceed with the new matrix, as 
described in [12]. The Neville elimination can be done without row exchanges 
if all the pivots are nonzero, as it will happen in our situation. The pivots pi^i 
are called diagonal pivots. If all the pivots are nonzero, then i = i 
and, by Lemma 2.6 of |12] 



P. 



det A[z - j + 1, . . . ,z|l, . . . ,j] 



detA[z - j + 1, 



i-lll. 



-J-l] 



1 < j < i < n. 



(2.2) 



The element 

TT^i j = 1 < J < j < i < n (2-3) 

Pi-hj 

is called multiplier of the Neville elimination of A. The matrix U := An is 
upper triangular and has the diagonal pivots in its main diagonal. 

The complete Neville elimination of a matrix A consists on performing the 
Neville elimination of A for obtaining U and then continue with the Neville 
elimination of U'^ . The pivot (respectively, multiplier) {i,j) of the complete 
Neville elimination of A is the pivot (respectively, multiplier) {j, i) of the 
Neville elimination of , if j > i. When no row exchanges are needed in the 
Neville elimination of A and U'^ , we say that the complete Neville elimination 
of A can be done without row and column exchanges, and in this case the 
multipliers of the complete Neville elimination of A are the multipliers of the 
Neville elimination of A if i > j and the multipliers of the Neville elimination 
of if j > i. 

A matrix is called strictly totally positive) if all its minors are positive. The 
Neville elimination characterizes the strictly totally positive matrices as follows 

my- 

Theorem 2.1. A matrix is strictly totally positive if and only if its complete 
Neville elimination can be performed without row and column exchanges. 
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the multipliers of the Neville elimination of A and are positive, and the 
diagonal pivots of the Neville elimination of A are positive. 

As it can be seen in [7] , the Said-Ball- Vandermonde matrices are strictly totally 
positive when the real numbers satisfy < ti < ^2 < • • • < tn+i < 1, and this 
fact has inspired our search for a fast algorithm, but this result will also be 
shown to be a consequence of our Theorem 3.2. 



3 Bidiagonal decomposition 



The Said-Ball basis iS„ = {sg (t), . . . , s^(t)} of the space n„(t) of the 
polynomials of degree less than or equal to n on the interval [0, 1] is defined 
by: 

sl^{t) = (L"/2J+^^i.(i _^)K2j+i^ < z < [(n- 1)/2J, 

S^{t) = (K2J+"-i)tK2j+l(l _ ^)n-i^ ^^/2J + 1 < z < n, 

and, if n is even 

where [m\ is the greatest integer less than or equal to m. 

From now on, we will call Said-Ball-Vandermonde matrices (SB- Vandermonde 
matrices in the sequel) the generalization of the Vandermonde matrices ob- 
tained when considering the Said-Ball basis instead of the power basis. The 
SB- Vandermonde matrices are therefore 



, Nli+i 

t2) 2 



/ / "-I \ , n + 1 / n-1 \ 

' (t)(i-^i)^ (t)(i 



n-l \ , n + 1 \ 



/ 'i-l I l\ , n + 1 

( — + l)t„ + l(l-t„+l) — 



/ \ n — 1 I , , 71—1 I -, , V 71 — 1 



n + 1 



(i-ti 
i\ "+i 

2 



'n-1 , -|\ " + 1 



(n-1 \ 11 



l\ "+i 
2 

2 



(t)« 



n-l\ — ^ 
2 ^ 



n + 1 
2 



n-1 

2 



n-1 , -ix n+l 

+ 1 1+ 2 



1 \ "+l 



(t) 
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in the case of odd n, and 



i-h) 2 



l-to 



n+2 



+2 
2 



(;)if(i-ti)^ (;)ti(i-t2)^ 



+2 

2 



«i' (l-«i) 

n + 2 



, nil. 



/n\ n+2 



n+2 
2 



( 



n-1 
2 



n— 2 n+2 

2 )in+l(l ~ tn+l)~ 



(i) 



l(1 -t„+i)2 



n+2 



n + 2 



1 )tn+l{^ " tn+l) 

n+2 

Qj^n+l 



in the case of even n. 



It must be observed that the SB-Vandermonde matrix A is the coefficient ma- 
trix associated with the following interpolation problem in the Said-Ball basis 
Sn'. given the interpolation nodes {ti : i = 1, . . . , 1} and the interpolation 
data {hi : i — 1, . . . ,n -\- 1} find the polynomial 

n 

Pit) = E (^ksiit) 

k=0 

such that p{ti) — bi ioT i — 1, . . . ,n + 1. 

From now on, we will assume < ti < ^2 < • • • < ^n+i < 1- 



Proposition 3.1. The determinant of the SB-Vandermonde matrix A defined 
above is 



det A = 
ifn is odd, and 

det^l 
ifn is even. 



n-l \ / w+l \ / Tt+3 ^ 



n -2\ n- V 
0/\l/\2 / '\^/\^, 



1 2 



n it,-ti), 

l<i<j<n+l 



'n\ / n+2 \ 

1 2 



'n-l' 

n-2 
V 2 / 



n 



n (ij-ti)^ 

2/ l<i<i<n+l 



Proof. Here we include the proof for the case in which n is odd. The proof in 
the even case is completely analogous. 
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Looking at [19], it can be observed that the matrix of change of basis from 
the Bernstein basis to the Said-Ball basis is a block-diagonal matrix M with 
triangular diagonal blocks, and whose determinant is 



detM 



(T)(T)(T)-te)te 



(«)(")■■■(::) 



(3,1) 



As it can be seen, for example, in |23] , the matrix of change of basis from the 
Bernstein basis 



Bn = {ht\t) 



n 



'l-tr-r, t = 0,...,n} 



to the power basis {1, t, t^, . . . , t"} is a lower triangular matrix of order n + 1 
whose determinant is 



det 



n\ In 



(3.2) 



Taking this into account, the matrix of change of basis from the power basis 
to the Said-Ball basis is MN~'^, and consequently, 

, , detM , 
det A = , , „ det V, 



detN 



where V is the Vandermonde matrix 



V 



1 ti ti 

1 t2 tl 



j-n \ 



1 i J.2 xn , 

Using the well-known formula for the determinant of a Vandermonde matrix 

l<i<j<n+l 

and the equations (3.1) and (3.2), the proof is concluded. □ 

The following two theorems will be essential in the construction of our algo- 
rithm for computing BT>{A) of a SB- Vandermonde matrix. 



Theorem 3.2. Let A = ( 

0'i,j)i<i,j<n+i a SB— Vandermonde matrix whose 
nodes satisfy < ti < ^2 < • • • < < ^n+i < 1- Then admits a 
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factorization in the form 



A — G1G2 ■ ■ ■ GnD ^FnFn-l ■ ■ ■ Fi 



-1 



(3.3) 



where Gi are upper triangular bidiagonal matrices, Fi are lower triangular 
bidiagonal matrices {i = 1, . . . , n), and D is a diagonal matrix. 

Proof. The matrix A is a strictly totally positive matrix (see [6lf7j ) and there- 
fore, by Theorem 2.1, the complete Neville elimination of A can be performed 
without row and column exchanges providing the following factorization of 
A-^ (see [HIS]): 

A"-^ = G1G2 ■ ■ ■ GnD^^FnFn-l ■ ■ ■ Fi, 

where F^ {I < i < n) are bidiagonal matrices of the form 



1 

1 



1 

-mi+i,i 1 

-mi+2,i 1 



V 



Gj (1 < < "ra) are bidiagonal matrices of the form 



GJ 



1 

1 



1 



(3.4) 



(3.5) 



y — m„+i^j 1 J 

and D is the diagonal matrix whose ith (1 < i < n + 1) diagonal entry is the 
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diagonal pivot pi^i — afj of the Neville elimination of A: 

D = diag{pi,i,P2,2, • • • ,Pn+l,n+l}- 



(3.6) 



First we obtain the expressions for the multipliers rriij and rriij, and for the 
diagonal pivots pi^i in the case of odd n. 

Taking into account that the minors of A with j initial consecutive columns 
and j consecutive rows starting with row i are 



n— 1 \ / n — 1 



detA[i,...,i + j-l\l,...,j]^^^j^ ^ j ^ 
(1-ti) 2 (l-tj+i) 2 • • • (1 - tj+j_i) 2 ni<fe<z<i+i-i(^; -4), 



if 3 < H^, and 



n-j 



det ^[^, . . . , 2 + J - 1| 1, . . . , j] = (X) f^^') ■ ■ ■ ( 

I\i<k<l<i+j-l{tl ~ tk) 



\n-j+l 



if j > ^ 



2 ' 



a result that follows from the properties of the determinants and Proposition 
3.1, and that rriij are the multipliers of the Neville elimination of A, we obtain 
that 



m. 



n+l 



(i-*»)^nLi(*'-*'-fc) 



J = 1,...,^; i^j + l,...,n + l, 



(l-tir-^+Hl-U-j)llL\iii-ii-k) ■ _ n+3 



(3.7) 

As for the minors of A^ with j initial consecutive columns and j consecutive 
rows starting with row i, they are: 



detA-[v..,i + i-l|l,...,i] = (X^-)(X ; V ...-2 

(1 - ti)^(i - h)"^ ■■■{!- t,)^fr'fif' ■ ■ ■ Ui<k<i<Ati - h), 



ifi<^andi + j-l<^, 
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if i < ^ and i + J - 1 > and 



detA'^[i,...,i + j-l\l,...,j] = ( 

n + l n+1 n-\-l 



Il^+n-i-j+2\ 
n— i+1 J \ n~i J \ n—i—j+2 J 



if i > 



These expressions also follow from the properties of the determinants and 
Proposition 3.1. Since the entries ruij are the multipliers of the Neville elimi- 
nation of , using the previous expressions for the minors of with initial 
consecutive columns and consecutive rows, it is obtained that 



/■ n—l 

2 ; 



rriij = < 



-^f ■ 7 = 2 1 = 1 7 — 1 



n+3. 



nLi(i-*^) 

n-t+2 1 

n — <+2 
^+n-i+2 1-tj 



; J = I,---, 



n+l 
2 ' 



2 ^,...,71 + 1; J — 1, . . . ,i ^, 



(3.8) 



Finally, the diagonal entries of D are: 



Pi,i 



n+3 



(3.9) 



The formulas for are obtained by using the expressions for the minors of 
A with initial consecutive columns and initial consecutive rows. 

As for the case in which n is even, proceeding analogously as in the odd case 
we obtain the following expressions for the multipliers rriij and ruij, and the 
diagonal pivots Pi^f 



m. 



(i-t,)^nLi(*'-*<-fc) 



J = l,...,f; 7 = j + l,...,n + l. 



(l-f,)"-i+l(l-ti-,)nLi(*'-*'-fc) _ n+2 



L (l-t»-l)"-J+2nL2fe-l-«i-/c) 



j = ^,...,n; 7 = j + + 

(3.10) 
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m 



i-i ''J' 

nLi(i-*'=) 

f +n-j+2 



2 'H- j = \ ?■ — 1 



„■ _ n+2. „■ _ 1 H 



n+6 
2 ' 



.'^ + 1; j = 1, ■ ■ ■ ,^ - 



ra+4 
2 ' 



(3.11) 



^,...,n+l;j = z-i^,...,^-l 



n+2 
2 



and 



-^j(i-t,)'^nfc<.(t.-4), ^ = i,...,f, 

I V n-i+l ) Y{lZ\{l-tu) ' 



(3.12) 



n+2 



2 ,...,n + l. □ 



Moreover, by using the same arguments of [21], it can be seen that this factor- 
ization is unique among factorizations of this type, that is to say, factorizations 
in which the matrices involved have the properties shown by formulae (3.4)- 
(3.6). 



Let us observe that the formulae obtained in the proof of Theorem 3.2 for the 
minors of A with j initial consecutive columns and j consecutive rows, and for 
the minors of with j initial consecutive columns and j consecutive rows 
show that they are not zero, and so the complete Neville elimination of A can 
be performed without row and column exchanges. Looking at equations (3.7)- 
(3.12) it is easily seen that rriij, rhij and pi^i are positive. Therefore, taking 
into account Theorem 2.1, this confirms that the matrix A is strictly totally 
positive. 

Theorem 3.3. Let A = ( 

'2jj)i<jj<n+i be a SB— Vandermonde matrix whose 
nodes satisfy < ti < t2 < • • • < ^n < in+i < 1- Then A admits a factorization 
in the form 

A = FnFn-l ■ ■ ■ FiDGi ■ ■ ■ Gn-lGn 

where Fi are lower triangular bidiagonal matrices, Gi are upper triangular 
{i = 1, . . . ,n), and D is a diagonal matrix. 

Proof. The matrix A is a strictly totally matrix [7] and therefore, , by Theorem 
2.1, the complete Neville elimination of A can be performed without row and 
column exchanges providing the following factorization of A (see [H] ) : 

A = FnFn-l ■ ■ ■ FiDGi ■ ■ ■ Gn-lGn, 
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where Fj {1 < i < n) are bidiagonal matrices of the form 



1 

1 



1 

mi+2,2 1 

Y 'm"n,n—i 1 y 

(1 < i < n) are bidiagonal matrices of the form 



GJ 



1 

1 



1 



(3.9) 



(3.10) 



V 

and D is the diagonal matrix 

D = diag{pi,i,P2,2, ■ ■ . ,Pn+l,n+l}- 

The expressions of the multipliers rriij (1 < i, j < n + 1) of the Neville elimi- 
nation of A, the multipliers rhij (1 < i,i < n + 1) of the Neville elimination 
of A^, and the diagonal pivots pi^i (l<i,<n-|-l)of the Neville ehmination 
of A are also in this case the ones given by Eq. (3.7) and Eq. (3.10), Eq. (3.8) 
and Eq. (3.11), and Eq. (3.9) and Eq. (3.12), respectively. □ 

It must be observed that the matrices Fj and Gi {i — 1, . . . ,n) that appear 
in the bidiagonal factorization of A are not the same bidiagonal matrices that 
appear in the bidiagonal factorization of A^^ , nor their inverses (see Theorem 
3.2 and Theorem 3.3). The multipliers of the Neville elimination of A and A^ 
give us the bidiagonal factorization of A and A~^, but obtaining the bidiagonal 
factorization of A from the bidiagonal factorization of A~^ (or vice versa) is 
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not straightforward. The structure of the bidiagonal matrices that appear in 
both factorizations is not preserved by the inversion, that is, in general, 
and (1 < i,j < n) are not bidiagonal matrices. See [14j for a more detailed 
explanation. 



4 The algorithm 



In this section we present a fast and accurate algorithm for computing BT>{A) 
for a totally positive SB-Vandermonde matrix A. Let us point out here that 
given A the matrix BV{A) represents both the bidiagonal decomposition of A, 
and that of its inverse A~^ (see Theorem 3.2 and Theorem 3.3). 

The algorithm will compute the multipliers of the Neville elimination of 
A, the multipliers rriij of the Neville elimination of A"^ and the diagonal pivots 
Pa of the Neville elimination of A, which are the entries of the matrix BV{A). 

We include here the algorithm for the case in which n is an odd number, the 
algorithm for the even case being analogous. 

The algorithm for computing the rriij given by Eq. (3.7) is: 

for i = 2 : n + 1 

n + l 

_ (l_t.)-2- 

^i.i — s+r 

(i-i,_i)^ 

for j = 1 : mm{i - 2, ^) 

end 
end 

for i = ^ -.n + l 

(1 — n + 3 " + 1 ) 

for J = : i - 2 



end 
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end 

The algorithm for the computation of the rfiij given by Eq. (3.8) is: 
f or i = 2 : ^ 



i-1 

for j — 1 : i — 1 



rriij = aux ■ tj 



end 

end 

mn + 3 -I = 

f or j = 1 : 



^-J + i tj(l-tj+l) ^-J 



end 



f or i = : n + 1 



f or j = 1 : i - ^ 



mi 



rriij = aux ■ int 



end 



f or j = i - ^ : i - 1 



int 



rriij — aux • int 
end 
end 

The algorithm for computing the diagonal pivots pi^i given by Eq. (3.9) is: 
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f or i = 1 : ^ 




aux — 1 
for k — 1 : i 

aux — {ti+i — tk) ■ aux 
end 

Pi+i,i+i = q-{l- tj+i)^ • aux 
end 

aux — 1 

f or A; = 1 : ^ 

aux = (1 — tk) ■ aux 
end 



_ _2_ 
aux 



aux 



= 1 



f or = 1 : 

awx = — tk) ■ aux 



end 



V n+S n+3 — 
2 ' 2 



• (1 — in+s) 



for i = ^ : n 



q = 



li^+n-i+l ' l-ti 



aux 



= 1 



for A; = 1 : i 



aux — {ti+i — tk) ■ aux 
end 

Pi+i,i+i = g • (1 - • aux 

end 

Looking at this algorithm is enough to conclude that: 

- The computational complexity of the computation of rriij, niij and pii, i.e. 
of the computation of BT>{A) is O(n^). 

- The algorithm has high relative accuracy because it only involves arithmetic 
operations that avoid inaccurate cancellation. 

- The algorithm does not construct the SB-Vandermonde matrix, it only 
works with the nodes {ti}i<i<n+i- 

As for the even case, the properties of the algorithm are exactly the same. 



5 Accurate computations with SB— Vandermonde matrices 

In this section algorithms for solving linear systems and for eigenvalue com- 
putation are presented for the case of a totally positive SB- Vandermonde 
matrix A. The algorithms are both accurate and efficient and are based on 
the algorithm presented in Section 4 for computing BT>[A). 

Let us observe here that, of course, one could try to solve these problems 
by using standard algorithms. However the solution provided by them will 
generally be less accurate since SB- Vandermonde matrices are ill conditioned 
(see the numerical experiments in Section 6) and these algorithms can suffer 
from inaccurate cancellation, since they do not take into account the structure 
of the matrix, which is crucial in our approach. 

5.1 Linear system solving 

Let Ax = 6 be a linear system whose coefficient matrix A is a SB- Vandermonde 
matrix of order n + 1 generated by the nodes {ti}i<i<n+i, where < ti < . . . < 
tn+i < 1- 

The following algorithm solves Ax — b in a fast way. 
INPUT: The nodes {ti}i<i<n+i and the data vector b e R"+^. 
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OUTPUT: The solution vector x G R"+^ 

- Stej) 1: Computation of BT>{A) by using the algorithm introduced in Section 
4. 

- Step 2: Computation of 

X = A-^h = G^G2 ■ ■ ■ GnD-^FnFn-i ■ ■ ■ F^b. 

Step 2 can be carried out by using the algorithm TNSolve of P. Koev |2Uj . 
Given the bidiagonal factorization of the matrix A, TNSolve solves Ax = 6 by 
computing the above matrix product. 

Although BT>{A) is computed with high relative accuracy, the accuracy of the 
solution vector will generally depend on the data vector b [23]. 

Taking into account that, as we have shown in Section 4, the computational 
cost of Step 1 is of O(n^) arithmetic operations, and the cost of computing 
whole product in Step 2 (from right to left) is also of 0(r2^) arithmetic oper- 
ations, the computational complexity of the algorithm for solving Ax = b is 
Oin"). 



5.2 Eigenvalue computation 

Let A be a SB-Vandermonde matrix of order n + 1 generated by the nodes 
{ti}i<i<n+ii where < ti < . . . < t„+i < 1. The following algorithm computes 
accurately the eigenvalues of A. 

INPUT: The nodes {ti}i<i<ri,+i. 

OUTPUT: A vector x G R"^^ containing the eigenvalues of A. 

- Step 1: Computation of BV[A) by using the algorithm introduced in Section 
4. 

- Step 2: Given the result of Step 1, computation of the eigenvalues of A by 
using the algorithm TNEigenvalues. 

TNEigenvalues is an algorithm of P. Koev [21] which computes accurate eigen- 
values of a totally positive matrix starting from its bidiagonal factorization. 
The computational cost of TNEigenvalues is of 0(?2^) arithmetic operations 
(see [2^) and its implementation in Matlab can be taken from |20]. In this 
way, as the computational cost of Step 1 is of O(n^) arithmetic operations, 
the cost of the whole algorithm is of O(n^) arithmetic operations. 
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6 Numerical experiments 



In this section we present two numerical experiments illustrating the accuracy 
of the two algorithms we have introduced in the previous section. 

Example 6.1. Let S15 be the Said-Ball basis of the space of polynomials with 
degree less than or equal to 15 in [0, 1], and let A be the SB-Vandermonde 
matrix of order 16 generated by the following nodes: 

1 1 2 317247 17 15 9 7 85 20 
— < — < — < — <-< — <-<-< — < — < — < — < — < — <-< — . 
16 13 11 13 4 18 5 9 15 30 26 13 10 11 6 21 

The condition number of A is: k,2{A) = 3.2e + 08. Let us consider the data 
vector 

b = (12, -3, 0, 1, 5, -7, 0, 2, 21, -4, 0, 9, -11, 6, -8, 0)^. 

We compute the exact solution Xe of the linear system Ax = b hj using the 
command linsolve of Maple 10 and we use it for comparing the accuracy of 
the results obtained in Matlab by means of: 

(1) The algorithm presented in Section 5.1. We will call it MM. 

(2) The algorithm TNBD of Plamen Koev [20j that computes BV{A) without 
taking into account the structure of A. 

(3) The command A\b of Matlab. 

In (2), the second stage in the solution of the linear system is the computation 
of the fast product (from right to left) of the bidiagonal matrices and the vector 
b. It is done in Matlab by using the same command as in (1): TNSolve of 
Koev [20]. 

We compute the relative error of a solution x of the linear system Ax = b hj 
means of the formula: 

II ^ — 1 1 2 

err = — . 

II h 

The relative errors of the solutions of Ax = b computed by means of the 
approaches (1), (2) and (3) are reported in Table 1. 



MM 


TNBD 


A\b 


5.1e-16 


2.2e-09 


3.9e-10 



Table 1 

Relative errors in Example 6.1 

Example 6.2. Let A be the SB-Vandermonde matrix of order 16 considered 
in Example 6.1. In Table 2 we present the eigenvalues Aj of A and the relative 
errors obtained when computing them by means of: 

(1) The algorithm presented in Section 5.2. We will call it MM. 
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(2) The algorithm TNBD |20] that computes BV{A) without taking into ac- 
count the structure of A. 

(3) The command eig from Matlab. 

In (2), the second stage in the computation of the eigenvalues is done in 
Matlab by using the same command as in (1): TNEigenvalues of P. Koev 
I20j. 

The relative error of each computed eigenvalue is obtained by using the eigen- 
values computed in Maple 10 with 50-digit arithmetic. 





MM 


TNDB 




1 Oe -1- 00 


4.4e - 16 


l.Oe - 12 


1.8e - 15 


9.4e - 01 


1.3e - 15 


2.1e - 11 


1.2e - 15 


7.0e - 01 


9.6e - 16 


2.5e - 11 


6.4e - 16 


5.2e - 01 


6.3e — 16 


1.3e - 11 


2.1e - 16 


3.1e - 01 


5.4e - 16 


7.9e - 12 


2.7e - 15 


1.4e -01 


1.3e - 15 


1.5e - 11 


1.3e - 15 


6.0e - 02 


5.7e - 16 


l.le - 11 


l.le - 15 


3.0e - 02 


4.6e - 16 


6.2e - 12 


4.6e - 16 


8.6e - 03 


4.1e - 16 


4.6e - 12 


1.3e - 14 


2.6e - 03 


9.9e - 16 


l.Oe - 11 


3.7e - 14 


6.1e -04 


5.4e - 16 


2.3e - 11 


7.2e - 14 


6.2e - 05 





l.Oe - 11 


3.1e - 13 


8.3e - 06 


4.1e - 16 


1.8e - 11 


6.4e - 13 


9.1e-07 


1.2e - 16 


4.6e - 11 


3.2e - 12 


5.5e - 08 


2.0e - 15 


1.2e - 10 


3.1e - 10 


5.0e - 09 


3.0e - 15 


2.3e - 09 


2.0e - 09 



Table 2 

Relative errors in Example 6.2 



The results appearing in Table 1 and Table 2 illustrate the good behaviour 
of our approach. In particular, the very different results obtained for the ap- 
proaches (1) and (2) show the importance of computing B'D{A) with high 
relative accuracy, since in both approaches the second stage is exactly the 
same. 
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For this specific matrix A the relative error obtained when computing the ma- 
trix BV{A) by using the algorithm we have presented in Section 4 is 2.8e— 15, 
while the relative error obtained when computing it by means of the command 
TNBD is 6.8e — 10. These relative errors have been computed for each solution 
B by using 

W B — Be II2 

err = — —— , 

II II2 

where B^. is the exact BV{A) computed in Maple 10. 
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